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ABSTRACT 

We give a new coherent description of the first-order Fermi acceleration of particles in 
shock waves from the point of view of stochastic process of the individual particles, 
under the test particle approximation. The time development of the particle distribu- 
tion function can be dealt with by this description, especially for relativistic shocks. 
We formulate the acceleration process of a particle as a two-dimensional Markov pro- 
cess in a logarithmic momentum-time space, and relate the solution of the Markov 
process with the particle distribution function at the shock front, for both steady and 
time-dependent case. For the case where the probability density function of the en- 
ergy gain and cycle-time at each shock crossing of the particles obeys a scaling law in 
momentum, which is usually assumed in the literature, it is confirmed in more general 
form that the energy distribution of particles has the power-law feature in steady state. 
The equation to determine the exact power-law index which is applicable for any shock 
speed is derived and it is shown that the power-law index, in general, depends on the 
shape of the probability density function of the energy gain at each shock crossing; in 
particular for relativistic shocks, the dispersion of the energy gain can influence the 
power-law index. It is also shown that the time-dependent solution has a self-similarity 
for the same case. 

Key words: acceleration of particles - shock waves - methods: analytical - methods: 
numerical - cosmic rays. 
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1 INTRODUCTION 

The first-order Fermi acceleration in shock waves is a widely 
known mechanism which generates non-thermal energetic 
particles in space. One of the most notable features of this 
mechanism is the power-law energy distribution of accel- 
erated particles in steady state. The mechanism works in 
various shock waves ranging from the earth's bow-shock to 
ultrarelativistic shocks associated with gamma-ray bursts. 
Near by the earth, satellites and space-crafts have directly 
observed energetic particles accelerated by this mechanism. 
Recent X-ray observations also discovered synchrotron X- 
rays from energet ic electrons in severa l supernova remnants 
(e.g. SN1006; see lKovama et al. Ill995t >. These electrons are 
considered to be accelerated in shock waves by this mecha- 
nism and have a power-law distribution. These observations 
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are believed to be evidence that cosmic-rays with energies 
below 10 15 eV (namely, 'knee') originate from supernova 
remnants in our galaxy. 

Since the basic theory was proposed in the late 1970's, 
this mechanism has b een investigat e d by numerous au- 
thors (for review see iDrurvl llQSfl iBlandford fc EichleJ 
Il987fl . While some of recent theoret ical interests are fo- 
cused on the non-linear problems JPrurv fc Voljd Il98ll : 
lEllison. Baring fc Jones] Il996t iBerezhko fc Ellison! I19991) . 
some of the linear problems, under the test particle ap- 
proximation, still remain to be clarified and need more 
close examinations, especially on the time development of 
the particle distribution and the acceleration in relativis- 
tic shocks. For non-relativistic shocks, this mechanism is 
described well by the diffusion-convection equation. Solv- 
ing this eq uation for steady state, the power-law solution 
is derived jAxford. Leer fc Skadronl Il977t iKrvmskvl Il977l : 
IBlandford fc Ostrikerlll978r) . Time-dependent solutions can 
also be examined based on this equation, usually with the 
aid of the Laplace transformation. However, to invert the 
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transform in analytical form is generally difficult and an- 
alytical solutions were derived only for several mathe mati- 
cally simple cases (e.g. iToptyghirJIlflSfl L iDrurvl l)l99lft pro- 
posed an approximation for more general situations by re- 
normalizing one of the analytical solutions using the first two 
cumulants of the p article distribution, w hich can be obtained 
in analytical form. iFritz fc Webbl |f 990) investigated the ef- 
fects of the synchrotron losses on the time-dependent so- 
lution for momentum-independent diffusion coefficients and 
calculated solutions by numerical inversion of the Laplace 
transforms. 

Recently, the particle acceleration in relativistic shocks 
has attracted some attention in relation to AGN jets, GRBs 
or ultra-high energy cosmic rays. In this situation, the prob- 
lem becomes more difficult because the anisotropy in the 
particle distribution is not negligible and, as a result, the 
diffusion-convection equation is no longer valid . Retur ning 
to the Boltzmann equation, iKirk fc Schneider! il987aT) de- 
rived the steady-state solution by means of a semi-analytical 
approach, performing the eigenfunction expansion. This ap- 
proach was rece ntly extended to ultrarelativistic shocks 
iKirk et al]l2000D . However, the time development still has 
not been treated in this way. 

There is an alternative approach established by iBelll 
( 1978) which is based on the acceleration process of individ- 
ual particles. From this approach, the acceleration process 
is described as a stochastic process of individual particles; 
particles are accelerated whenever they repeat the cycle of 
crossing and re-crossing of the shock front, where the en- 
ergy gain per one-cycle and the cycle-time are both regarded 
as stochastic variables. The description needs not introduce 
the assumption of the isotropy on the particle distribution. 
It can therefore be n aturally exten ded to the acceleration 
in relativistic shocks <|Peacocklll 98lT) . In our previous paper 
iKato fc Takaharall200fir we investigated the acceleration 
process in relativistic shocks from this approach utilizing 
random walk theory. However, because the previous studies 
dealt only with the energy gain per one-cycle and not with 
the cycle-time, they are restricted within the steady prob- 
lems and can not treat the time development of the distri- 
bution of particles, even for non-relativistic shocks. Further- 
more, the derivation of the distribution function of particles 
from this approach and that of the power-law index are still 
not satisfactorily established. 

Monte Carlo simulations, the third approach, can 
make a direct estimation of the distribution func- 
tion even for relativist ic and ultrarelativi s tic shocks 
jKirk fc Schneider! Il987bl jEllison et alj Il99(j lOstrowskil 
Il99ll: iBednarz fc Ostrowskillf996llf998f) . However, it needs 
a large-scale simulation to obtain sufficiently accurate re- 
sults, and physical interpretations of the results in terms of 
analytical models are also desirable. 

In this paper, we reanalyse the acceleration process 
based on a stochastic, single-particle approach. By treating 
the cycle-time as a stochastic variable explicitly, our method 
can describe the time development of the particle distribu- 
tion function. The method is also applicable to relativistic 
and ultrarelativistic shocks as well as non-relativistic shocks. 
The paper is organized as follows. In Section|5| we formulate 
the mechanism as a stochastic process. In Section |3] for a 
case where the property of the one-cycle of the shock cross- 
ings has a scaling law in momentum, we investigate steady 



and time-dependent solutions of the distribution function of 
particles by an analytical way. In Section 2] as a check, we 
apply this theory to non-relativistic shocks and compare the 
results with conventional results. In Section [K] we consider 
relations between the one-cycle properties and the result- 
ing particle distribution function for both steady and time- 
dependent case. Conclusions are given in Section |S| 



2 PROBABILISTIC DESCRIPTION OF THE 
ACCELERATION PROCESS 

In this section, we formulate the acceleration process of a 
particle as a stochastic process, and then relate it with the 
distribution function of particles. We also give numerical 
methods for obtaining the solution. 

2.1 Basic acceleration process in shock waves 

From the single particle point of view, the basic mechanism 
of the particle acceleration is explained as follows. First, be- 
cause of scattering caused by magnetic irregularity existing 
in the plasma in both sides of the shock front, particles move 
like random walk and, as a result, a fraction of the particles 
can repeat crossing and re-crossing of the shock front many 
times; this can be modelled by that in each cycle of the shock 
crossing, a particle in the downstream region returns to the 
shock front at a returning probability P rc t, and otherwise 
escapes from the acceleration region. Second, since the elec- 
tric field approximately vanishes in the respective plasma (or 
'fluid') rest frames in both regions, the energy of a particle 
measured in the respective rest frames is unchanged while 
the particle stays in one of the regions. Consequently, in a 
fixed reference frame, particles gain energy whenever they 
repeat the crossing cycle, owing to the difference between the 
fluid speeds; in particular for sufficiently relativistic parti- 
cles, letting the relative fluid speed between the upstream 
and downstream region be V le \(> 0) and its Lorentz fac- 
tor be r rc i = (1 — V^-ei/c 2 ) -1 / 2 , we can calculate the energy 
gain of a particle for one-cycle of the shock crossings, that 
is, a downstream to upstream to downstream cycle, by per- 
forming two successive Lorentz transforms together with the 
above condition: 

pt = r? cl (l - Yf cosflr) (l + cos#) p u (I) 

where pi and pt are the magnitudes of momentum of the 
particle before and after the cycle, respectively. 6\ and 8' 2 
are the angles at the shock crossing from the downstream 
to upstream region measured in the downstream rest frame 
and from the upstream to downstream region measured in 
the upstream rest frame, respectively. Considering the pos- 
sible ranges of these angles to cross the shock front, it is 
straightforwardly shown that the energy always increases at 
every shock crossing cycle (here, energy loss mechanisms are 
not taken into account). Because these angles are regarded 
as stochastic variables on account of the scattering process, 
the energy gain is also regarded as a stochastic variable. Ow- 
ing to the diffusive motion of the particles, the cycle-time is 
also regarded as a stochastic variable. 

This is a basic picture of the acceleration mechanism 
from the single particle point of view. Fig. illustrates the 
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Figure 1. A schematic illustration of the acceleration of a par- 
ticle in a non-relativistic, parallel shock wave. 



one-cycle of the acceleration process of a particle in a non- 
relativistic, parallel shock. 

From a theoretical viewpoint, the problem of the ac- 
celeration process consists of two parts. One is how parti- 
cles are accelerated for one-cycle of the shock crossings. The 
other is how the distribution function of the particles evolves 
when the one-cycle properties are given. The former is de- 
termined by details of the scattering process and is quite a 
difficult problem unless the diffusion approximation is appli- 
cable; the properties of both turbulent magnetic field in the 
vicinity of the shock front and transport of particles in such 
turbulent fields themselves remain long-standing problems. 
On the other hand, it is possible to establish a well-defined 
theory formally for the latter problem. In the present paper, 
we therefore construct that theory from the point of view of 
the single particle approach. 

2.2 Description as a Markov process 

In this paper, we investigate the acceleration process at the 
shock front; our primary aim is to solve the steady state 
and time-dependent solution of the distribution function of 
particles at the shock front. Therefore, we concentrate only 
on the state of particles, (p,t), at the crossing of the shock 
front, where p is the momentum of a particle and t is the 
time respectively; we measure p in the downstream frame 
and t in the shock rest frame, respectively. The history of 
acceleration process of a particle observed at the shock front 
can be described by a sequence of the state of the particle 
at each end of the acceleration cycle, that is, at each shock 
crossing from the upstream to downstream region: 



(Poi*o),(Pi,ti),(p 2 . i 2) 5 



(2) 



Such a sequence is terminated when the particle escapes 
from the acceleration region, usually from the downstream 
region. We consider this process in the p— t space where each 
state of the sequence defines a point, which is called 'state 
point' in the following. If the changes in p and t at each cycle, 
(Ap, At), are stochastic variables as already mentioned, the 
acceleration process can be regarded as a stochastic process 
in the p — t space. Usually, this is a Markov process because 



the changes (Ap, At) mostly depend only on the momentum 
at the last cycle. 

Although the history of (p, t) provides complete infor- 
mation of the acceleration process at the shock front, we 
consider only the history of (p, t) in this paper, where p is 
the magnitude of the momentum, because only p concerns 
most practical purposes. In order to treat the history of (p, t) 
solely as a Markov process, we also assume that the prob- 
ability density of the changes in p and t at one-cycle does 
not depend on the direction of the momentum at the last 
cycle, at least approximately. (If this assumption is not ap- 
propriate, we must deal with the state of a particle by (p, t) 
or (p,p, t), where fj, is the cosine of the angle between the 
shock normal and the momentum of the particle.) For con- 
sidering the acceleration of relativistic particles, it is better 
to describe this process in a logarithmic momentum space, 
because the changes in p at each cycle are usually multiplica- 
tive rather than additive as shown in equation Q. Further- 
more, since we consider only mono-energetic injection with 
a certain momentum of po measured in the downstream rest 
frame in the fundamental part of the following description, 
we use a new quantity 



q := ln(p/p ) 



(3) 



instead of p to simplify the description. 

Thus, the acceleration process of a particle can be con- 
sidered as a two-dimensional Markov process on the q — t 
plane. This stochastic process is described by a probability 
density function which describes the transition on the q — t 
plane for one cycle of the acceleration, p(Aq, At; q), where 
Ag and At are respectively the changes in q and t at one- 
cycle; q in this notation denotes the logarithmic momentum 
at the last cycle representing the Markov property of this 
process. (We also use the term 'step' instead of 'cycle' for 
the point of view of Markov process in the following.) Since 
a particle can be lost from the acceleration region at each 
cycle at a probability of P esc (q), this density is a defective 
one; the integration of it over all area is not normalized to 
unity, but to the return probability, P rc t(g) = 1 — Pesc(fj'), as 
follows: 



dAq / dAt p(Aq, At; q) = P ret (q). 



(4) 



In the present paper, for simplicity, we assume Aq > 0, that 
is, any energy loss mechanisms are not efficient. (However, 
to generalize the following description to include the possi- 
bility of Ag < is not difficult.) The correlation between the 
energy gain Ag and the cycle-time At, which can be con- 
sider able for relativistic shocks (see iBednarz fc Ostrowskil 
1996), can be included in the functional form of p. For later 
convenience, we define here the moments of Ag and At as 
follows: 



dAq / dAt (Ag)" p(Aq,At;q), (5) 



<(Ag)") 



<(Atr>:=— / dAq I dAt (At)" p(Aq,At;q), (6) 

"ret 







(I 



a Aq := y/({Aq)*) - (Ag)*, 
a At := V((At) 2 ) - (At)2. 



(7) 
(8) 
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These quantities are generally dependent on q at the last 
cycle. As already mentioned, although the functional form 
of p(Aq, At; q) is determined by microscopic physics, we do 
not specify it and discuss the process formally in this paper 
(we adopt only an approximate model of p in Section 2] and 
simple toy-models in Section 

Defining the probability density function of the state 
points at nth step, p n , the progress from (n — l)th step to 
nth step is expressed by 



Pn(q,t) 



dt' p(q-q',t-t';q')p n -i(q',t'). ( 9 ) 



Again, the integrals of these p n 's over all area are not nor- 
malized to unity, but to the probability at which a particle 
continues the acceleration cycle at least n times, P n : 



dq 



dt p n (q,t) = P n - 



(10) 



Since we consider now the case where particles are all in- 
jected at the same state of (go, to) = (0, 0), we have 



pi(q,t) = p(q,t;q = 0). 



(11) 



The integral in equation © includes a plain convolution 
with respect to t. It is better to represent such equations in 
the form of the Laplace transform with respect to t. Intro- 
ducing the following notation 



f{s) 



f(t)e- st dt, 



(12) 



where f(t) is a function of t, equation © is transformed to 



Pn(q, 



p(q 



,s;q)p n -i(q',s) 



(13) 



In a mathematical sense, the above stochastic process 
can be regarded as a terminating renewal process in the 
probability theory (cf. lFelleJl97ltlCox fc Milleill965lh How- 
ever, since the present process is two-dimensional and also 
has a Markov property, it is a more complex problem. Fur- 
thermore, the stochastic process defined by equation 1131 for 
a fixed s can be regarded a s the random walk proc ess dealt 
with in the previous paper dKato fc Takaharall200lh . except 
that the present process is in q-space with the Markov prop- 
erty, always increases the value of q and has no absorbing 
barriers. 



2.3 The density of state points 

A history of acceleration process of a particle can be repre- 
sented by plotting the state points on the q — t plane. When 
such plots are made for many particles and superposed as 
shown in Fig. [5] we can define the density of the state points: 



9(q,t) ~J2p n (q,t), 



(14) 



where the density is normalized by the number of the parti- 
cles. This is a key concept to treat the present problem. The 
physical meaning of the product ^{q,i)dqdt is the expecta- 
tion value of number of particles which are initially injected 
at (q, t) = (0, 0) and then cross the shock front from the 
upstream to downstream side within ranges [q, q + dq] and 




10" 



Figure 2. A scatter plot of the state points superposed for many 
particles obtained by Monte Carlo simulations (see Section lil for 
details). The unit of time is taken the mean cycle-time for q = 0. 
The density of these points defines ^f(q, t) in equation 1141 . 



[t, t+dt], per unit injection. We can therefore relate this den- 
sity with the flux and distribution function of the particles 
at the shock front (see the following subsection). 

Combining equations © and 11411 . we obtain the fol- 
lowing integral equation: 

9(q,t) = p(q,t;0) + f dq' f dt' p(q-q' ,t-t' '; </)*(?' ,*'). (15) 
Jo Jo 

This equation determines $ from only p, and is a fundamen- 
tal equation to our description. Taking the Laplace trans- 
form of this integral equation with respect to t, we obtain 

*(<?, s) = p(q, s; 0) + / p(q-q',s;q')$(q',s) dq'. (16) 



For a fixed s, this integral equation is regarded as a Volterra 
equation of the second kind. Since it is generally difficult to 
solve this equation in analytical form, we will give numerical 
methods in Section 12.51 In the following, we will use the 
following relation 



(17) 



where *o(?) := *(?,0). 

Some quantities related with the acceleration process 
are described in terms of VP as follows. From equations (1101 
and 114H . we obtain 



dq / dt ty(q,t) = 
Jo J n 



Mq) dq = ^P n = (n), (18) 



where (n) is the expectation value of the number of cycles in 
an acceleration process. (The probability at which a particle 
continues the cycle just n times before escaping is given by 
Pn — Pn+i-) The probability at which a particle escapes from 
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the acceleration region before attaining to q for < t < oo 
is given by 



Mq) 



P csc (0)+ / dq'P csc (q') / dt'y(q',t') 



Pcsc(0)+ / P csc (q')^ (q')dq' 

Jo 



(19) 



because a state point at q does not make the next state 
point at the probability P csc (g). We can show -Pe(?) -» 1 as 
q — > oo. The probability at which a particle attains to q for 
< t < oo is of course given by -Pa(?) = 1 — Pe((?), and we 
obtain from equation 119H 



dPA(q) 
dq 



Pesc (<?)*()(<?) 



(20) 



Note that -Pa(?) does not approach unity as q — > but to 
1 — Pcsc(0) because ^ includes no contribution from the state 
points of the particles just injected at (0, 0), by its definition. 
The mean acceleration time for attaining to q is defined by 



f~t*(q,t)dt 



ds 



/*o(«). 



(21) 



From the point of view of the renewal process, the den- 
sity of state points ^ and the integral equation 1151 cor- 
respond to the renewal density and the renewal equation, 
respectively. A con cept similar to ^ was al so introduced in 
the previous paper jKato fc Takaharal200ll) to deal with the 
random walk of particles. 



2.4 Particle distribution function at the shock 
front 

Here, we relate the density of state points introduced in the 
previous subsection with the particle distribution function at 
the shock front. We consider the acceleration of relativistic 
particles in a plane shock, where the velocity of the particles 
can be approximated by the speed of light c in the shock rest 
frame and in both fluid frames. 

Consider first an impulsive injection of unit particles 
with q = at t — at the shock front. Letting /j, be the cosine 
of the angle between the shock normal and the direction of 
particle momentum measured in the shock rest frame, we 
can define the distribution function of particles at the shock 
front for this injection, W s {[q, fi,t), where q is measured in 
the downstream rest frame; \x and t are measured in the 
shock rest frame. In terms of this function, the one-sided 
flux of particles with logarithmic momentum q at the shock 
front from the upstream to the downstream is represented 
by 

S Bf+ {q,t)= f c l xW si {q,iJ l ,t)d l x = c{n)+W s{+ {q,t), (22) 
Jo 

where we introduce the following notations: 



W s{+ {q,t):= / W st {q,n,t)d^, 



<M>- 



tiW s {(q,fi,t)dfM j Wsi+(q,t), 



(23) 



(24) 



where + is, in general, dependent on q . On the other 
hand, as already mentioned, because 



${q,t)dt = S st+ (q,t)dt (25) 
for the present situation, we obtain the following relation 



W s( +(q,t) = 



1 



(26) 



Because this W s t+{q, t) can be regarded as a Green's function 
of the time-dependent solution of the distribution function 
at the shock front, the solution for more general injection 
rate Q(t), F B { + (q,t), is constructed as 



F a( +{q,t) 



W si+ (q,t~t')Q(t')dt' 



V{q,t-t')Q(t')dt', 



(27) 



C W + Jo 

where F s f + is defined from F s f like W s i+, or in the form of 
the Laplace transform 

1 



F*f+ (q, s) 



c(/x). 



-*(q,s)Q(s) 



(28) 



In the following, we deal with only the steady injection, 
that is, Q(t) = Q = const for t > and Q(t) = for t < 0. 
For the case, the steady state solution (t — > oo) is given by 



F s(+ (q, oo) = 



-*o(g), 



(29) 



using the relation l|17^ . In the following, to simplify the nota- 
tions, we also express the time-dependent solution in terms 
of the 'cut-off function' defined by 



Q(q,t) 



J*9(q,t')dl/ _ J^jq^dt' 



r o °°*( g ,t')dt' * (q) 

resulting in 

F al+ (q,t) = F al+ (q,oo)e(q,t). 



(30) 



(31) 



The Laplace transform of this function with respect to t is 
written by 



Q(q,s) 



s *o(q) ' 



(32) 



Because ^(q,t) always takes non-negative value, we can see 
that < Q(q,t) < 1 for any values of (q, t) for the steady 
injection; the time-dependent solution always lies under the 
envelope of the steady state solu tion. This fun ction is iden- 
tical to the function <fi(t,x,p) of iDrurvl (| 199 if ) at the shock 
front (x — 0) [see equation (9) of his paper] , which was used 
in the analysis of time-dependent solutions of the diffusion- 
convection equation. 

It should be noted that, in some cases, the steady state 
solution _F a f+(q,oo) can be evaluated without introducing 
the density function Substituting equation 1201 into equa- 
tion 12911 . we obtain 



F a t+(q, oo) 



c(n) + P BSC (q) 



dP A (g) 
dq 



(33) 



Therefore, if -Pa(?) is estimated by another way, the steady 
state solution is obtained immediately; conventional single- 
particle approaches can estimate Pa (q) approximately (e.g. 
lBelil978t1PeacoclJl98lh and therefore can obtain the steady 
state solution. Although, to treat the time development of 
the distribution function, we must treat the function ^ ex- 
plicitly. 
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2.5 Numerical Methods 



and 



By numerical methods, the integral equation 11611 can be 
solve d straightforwardly (see chapter 18.2 of IPress et alJ 
2002). For a given range < q < g ma x, dividing the range 
into N meshes with uniform intervals of h := q m a*/N and 
adopting the trapezoidal rule to the integral in the equation, 
and introducing notations qt — ih, pij := p(qt — qj, s; qj) and 
tyi := ty{qi, s) for i = 0, 1, N (we omit the dependence on 
s for simplicity here), equation l|16[l is discretized as follows: 



Poo, 



i h ~ 



9i 



f i-i x 



(34) 



(35) 



where i = 1, 2, A r . Thus, ^i's are trivially solved by for- 
ward substitution. Note that the interval h must be fine 
enough to resolve the kernel p. 

Once ^(q, s) is obtained by this way for each s, we can 
then make a numerical inversion of Laplace transform de- 
scribed in Appendix lAlto obtain ^{q,t). We can also obtain 
F a {+(q, t) for general injection or Q{q,t) for steady injection 
numerically by inverting equation 1281 or 13211 . respectively. 



3 SELF-SIMILAR SOLUTION OF PARTICLE 
DISTRIBUTION 

One of the most well-known features of the first-order Fermi 
acceleration is the power-law distribution function of accel- 
erated particles in steady state. Theoretically, it is just a 
consequent of some scaling laws, or similarities, in the one- 
cycle probability density function p; in theoretical studies on 
the first-order Fermi acceleration, such a scaling law is usu- 
ally assumed implicitly or explicitly. The scaling-law may 
realize when, for example, the power spectrum of the turbu- 
lent magnetic field has a power-law feature. The fact that the 
power-law distribution is really observed in various shocks 
may indicate such a scaling law often realizes well. In the 
following, we also examine such a case in an analytical way. 



3.1 One-step probability density function with a 
scaling law 

We consider here the case in which the transition of each 
step (Aq, At) satisfies the following scaling law; for normal- 
ized variables Aq and At/(p/po) a (= e~ aq At), where a is 
a constant parameter, the probability distribution function 
defined for these variables is independent of q at the last step 
or, in other words, common for all q. This means that the 
one-step probability density p(Aq, At; q) can be represented 
by a single function <p(Aq, e~ aq At) as 



p(Aq,At;q)=e- a ^(Aq,e- a «At) 
and 

/5(A<?,s;g) = <p(Aq,e aq s) . 



(36) 



(37) 



The similarity of this function is clear. Note that P rct (g), 
Pcsc(q) and the moments of Aq are all independent of q: 



fret = / <po(Aq)dAq = const 
Jo 



(38) 



((Ag)"> 



Pre 



(Ag) n i^o(Ag)ciAg = const, 



(39) 



where tpo(Aq) := tp(Aq,s — 0). The integral equation ill HI 
is reduced to 

*(g, a) = <fi(q, s) + [ q>{q - q' , e aq 's)^(q', s) dq' . (40) 
Jo 

3.2 Steady state solution 

3.2.1 Power-law solution 

We first consider the steady state solution here. Setting s = 
in equation Q4O0 . we have 

*o(g) = Ml) + [ M<1 ~ q')^o(q') dq'. (41) 
Jo 

Again, denoting the Laplace transform with respect to q by 



f W 



f(q)e q dq, 



(42) 



where f(q) is a function of q, we can write the exact solution 
of equation 1411 in the form of the Laplace transform as 
follows 



(43) 



i-£5(0) - 

As well known, the inversion of this transform is given 
through the Bromwich integral on the complex S-plane: 



1 - <p* {6) 



(44) 



where 7 is a real constant taken so that the vertical line 
Re(f9) = 7 lies on the right of the all poles of (p ( *(f9)/[l — 
<Po(9)]. While the numerical solution can be obtained di- 
rectly (see Appendix [XJ, we investigate the asymptotic be- 
haviour of equation 1 IH here. 

Assuming that if 0(6) — > as \9\ — > 00 (we expect this 
condition usually holds), the contour of the integral in the 
last equation is closed through the left of the #-plane. The 
integrand obviously has poles at the points where the con- 
dition <f> 0(8) = 1 is satisfied. Because ifio(Q) = P rc t < 1 and 



Aqf (Aq)e- eAq dAq < 0, 



(45) 



the equation (p(>(6) — 1 has the unique negative root on the 
real axis of the complex (9-plane. We write this root as 6-. 
As q becomes sufficiently large, because only the pole at 9- 
dominantly contributes to the integral 1441 . we obtain the 
asymptotic solution (q — > 00) 

*o(?) ~ A^e-*, 



where A 



A:=- 



d8 



(A > 0) and 

Aqfo(Aq)e XAq dAq. 



(46) 



(47) 



The index A is, by definition of 9- , determined as the unique 
positive root of 



fo{Aq)e XAq dAq = 1. 



(48) 
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Recalling q — ln(p/po), the steady state solution 1291 for 
large q is written by 

4^,00)^^1 (JL) , (49) 

where F^(p)dp = F B {+(q — \n(p/po))dq. This is the power- 
law solution whose index is given by a — 1 + A. 

3. 2. 2 Approximations 

If the dispersion of Aq, tJAq, is negligible in equations 14811 
and 147H . namely (ta 9 <C A -1 , we can approximate <fio(q) = 
PvctS(q — (Ag)) there, resulting in 

In P rct 



A = — • 



(Ag) 



(50) 



This expression is equivalent to equation (5) of IPeacockl 
il98lf) . We will mention the validity of this formula again in 
Section [3] The condition to use this approximation can be 
rewritten as 



> - In P rc 



(51) 



0~Aq 

Furthermore, if <ta 9 <C (Aq), we have 

A = (Ag). (52) 

On the other hand, if (Ag) < A -1 and ua ? <C A -1 , expand- 
ing exp(Aq) to the first-order of Ag in both equations 14811 
and 1471 1 . we obtain 

A= p P ; B ; r , A = P rct <Ag). (53) 

As shown in Section 0] this approximation is applicable to 
the acceleration in non-relativistic shocks. 

3.3 Time-dependent solution 

3. 3. 1 Self-similarity 

Here, we show that the time-dependent solution of equation 
l|40|l has a self-similarity for q 3> (Ag). Firstly, for g 2> (Ag), 
we have approximately 



*(g,s) 



<f>(q - 9 e aq s)^(q',s)dq. 



(54) 



This equation must be regarded as a relation for large g, not 
as an integral equation defined over all region of g; such an 
integral equation would have only a trivial zero solution. For 
a positive constant a, it is seen that the function defined by 



* a (g,s) := *(g + lna,a a s) 



(55) 



also satisfies the above relation instead of in other word, 
the relation does not change its form under the transforma- 
tion q — > q + In a and s — > a~ a s simultaneously. Secondly, 
the steady state solution 1461 , which is regarded as a 'bound- 
ary condition' at s — 0, is only multiplied by the constant 
a~ A under this transformation. Finally, because the relation 
15411 is linear, we obtain the following self-similarity of the 
density of state points for q 3> (Ag): 

\i/(g + In a, a~ a s) = a~ x ^/(q, s), 



V(q + lna,a a t) = ffl - (a+A) *(g, t). 



(56) 
(57) 



For the cut-off function, we obtain 
B(g + In a, a a t) = Q(q,t). 



(58) 



Defining functions for p instead of g, \f , ' p -' (p)dp := *l/(g = 
]n(p/po))dq and 0^ p '(p) := 0(q = ln(p/po)), we obtain 



<H ip \ap,a a t) = a- (1+a+X) ¥ p) (p,t), 
Q M (ap,a a t) =9 (p) (p,t). 



(59) 
(60) 



For later convenience, we derive a relation between the 
steady state solution and the time-dependent solution here. 
Because ip(q, i) is usually a concentrated function of t, the 
Laplace transform of it, y>(g, s), has a characteristic value 
of s = s» so that (p{q,s) can be approximated by (po{q) for 
s < s*. Therefore, in the region of q — s plane where the 
condition e aq s < s* is satisfied, the integral equation 1401 is 
approximated by 

*(g, s) = ^o(q) + / Mq - q'Mq, (fii) 

Jo 

resulting in ^(g^) ~ ^o(g) in that region. Thus, we can 
define 



9.W 



ln(s,/s) 



(62) 



so that we have ^(g, s) ~ 'l'o(g) for s < s* and g < g„(s) 



3.3.2 Approximate solution 

When $ varies with g much more slowly than (p, we can 
approximate in the integrand of equation 14011 as 

*(g', s) ~ *(g, a) + (q 1 - 9)§^(g, »)• (63) 

A rough criterion to use this approximation may be given 
by (Ag) <C A -1 and ua, <C A -1 , since the typical variation 
scale of is given by A -1 for the steady state solution l|46|l . 
Thus, for q S> ( Ag) , the integral equation I40H is reduced to 



9* 

*(g, s) = Bi(q, s)*(g, s) - B 2 (q, s)-g-{q, s), 



where 

B 1 (q,s) := 
S 2 (g,s) := 



<f>(Aq, e aq s)dAq, 



Aq<p(Aq,e°"'s)dAq. 



(64) 

(65) 
(66) 



For a fixed s, equation l|64ft can be regarded as an ordinary 
differential equation for and the solution is given by 



*(g, s) = C(s) exp 



i-£i(g',s) 
S 2 (g',s) 



dg 



(67) 



where C(s) and go(s) are functions of s; we can choose the 
function go(s) arbitrarily here. If we choose go(s) = g*(s), 
which was defined in equation 1621 . we can determine C(s) = 
*o(g*(s)), because *(g,s) = *o(g) for g < g*(s). Thus, the 
solution is written by 



*(g,s) = * (g,(s))exp 



l-£i(g',s) 
B 2 (g',s) 



f/g' 



(68) 
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It is easily seen that the above expression satisfies the simi- 
larity 15611 . because Bi(q + In a, a~ a s) = Bi(q,s), B2(q + 
lna,a~~ a s) = B2(q,s) and q*(a~ a s) — q*(s) + Ina. For 
q < g*(s), because Bi(q,s) = P rct and B 2 (q,s) = P rct (Ag), 
we obtain 



9 1 



Pret{Aq) 



(g- g„(s)), 



(69) 



which is consistent with the steady state solution 1461 with 
A of equation 1531 . For q ^> (Ag), where the asymptotic 
solution 1461 is applicable for steady state, we finally obtain 
the following asymptotic expression 

*(g, s) = yrV r(!M) , 



(70) 



(71) 



where 

f q l-B 1 (q',s) J , 
I(q,S)~ / , , r dq. 

J B 2 {q',s) 

We will use this approximation in the following section. Al- 
though it is usually difficult to invert the above expressions 
analytically, it is able to do it numerically (see Appendix 
U\l. Under the above approximation, the mean acceleration 
time defined by 12 U is represented by 



-(At) , 



(72) 



aP let {Aq) 

where {At)o denotes the mean cycle-time (At) for q — 



4 APPLICATION TO NON-RELATIVISTIC 
SHOCKS 

In order to confirm that our method surely reproduces the 
well-known results, we apply the method to the acceleration 
in non-relativistic shocks in this section. We consider a sim- 
ple one-dimensional parallel shock as in Fig. with shock 
speed V s h(<C c). The fluid speeds in the upstream and down- 
stream region measured in the shock rest frame are denoted 
by V u and Vd, respectively, and are assumed uniform in each 
region. The compression ratio is given by r = V u /Vd- 



4.1 One-step probability density function 

When we apply our method to a specific case, we must first 
specify the one-step probability density p. Assuming that p 
obeys the scaling law 1361 discussed in the previous section, 
we derive here an approximate expression of y> applicable to 
non-relativistic shocks, by a practical way. 

First, because the changes in q and t at one-cycle are ex- 
pected to be mutually independent, the dependences of ip on 
Ag and At can be separable as <p(Aq, At) = ip q {Aq)ip t (At) . 
Since the particle distribution function at the shock front is 
approximately isotropic, we have 



(Ag) 



4K 
3 



V d 



O Ag 



IK 
3 



(73) 



c ' 6 c 

[cf. equation Q, or equation (7) of lBelili978| . These two 
quantities are expected to be small enough, compared with 
A -1 , to use the approximations described in Sections 13.2.21 
and 13.3.21 Therefore, using these approximations, we can 
proceed without specifying the functional form of tp q here. 
The time dependent part (p t can be approximated as 



follows with the aid of a solution of the diffusion-convection 
equation and results of Monte Carlo simulations. In a one- 
dimensional diffusion process with a diffusion coefficient D 
in a moving medium with speed V, the time distribution 
function for returning from a distant point at x = a (a > 0) 
to a barrier located at x = is given by 



9(f) 



V^nDt 3 



exp 



(a + Vtf 
ADt 



(74) 



or in the form of the Laplace transform 

3( S )=exp{-^[y + (y 2 + 4D S ) 1 / 2 ]} (75) 

fsee lCox fc MilleJll965l:lLagage fc Cesarskvlll983Fl . where V 
can be positive or negative and, for the present problem, V = 
— VIi for the upstream region and V = Vd for the downstream 
region. Since this function already includes the influence of 
the escape of particles, the return probability is given by 



P fet (a)=s(0) 



1 for V < 

exp(-aVVD) for V > 0. 



First two moments of the return time t are given by 



(t) 



(t 2 ) - (ty 



2aD 
W 



(76) 



(77) 



However, the above results are not applicable directly 
to calculate the cycle-time; in order to do this, we need 
the distribution of the residence time, that is, the time be- 
tween entering one of the fluid regions and leaving the region 
through the shock front, in each region. Here, we employ 
Monte Carlo simulations to obtain an approximate expres- 
sion for the residence time distribution. It is known that, 
in non-relativistic shocks, the properties of particle motion 
have little dependence on the details of the scattering pro- 
cess. Thus, we can utilize the large-angle scattering model, 
which was dealt with the previous paper jKato fc Takahara 
l200ll) . in the simulation to estimate the residence time for 
isotropic injection of particles at the shock front. Performing 
the Monte Carlo simulations, we found that the distribution 
of the residence time for isotropic injection at the boundary 
is approximated well by g(t) in equation 1741 if we set the 
parameters for the diffusion process as follows: 



a = 3 CT °' 



n 1 2 

D = TO, 



(78) 



where to is the mean free time of the particles, defined for 
the large-angle scattering model, measured in the fluid rest 
frame. Thus, we can use g(t) with a — 4D/c as an approx- 
imation of the residence time distribution for each region. 
Taking the convolution between the residence time distribu- 
tions for the upstream and downstream region, we finally ob- 
tain the distribution of the cycle-time in the form of Laplace 
transform 



ifit(s) ~ G(s) = exp 



{-2[^ d 



+ ("u + 



4D U 



,1/2 



/ 2 
(^d - 



4D d 



(79) 



where u u := V u /c and v& ■= Vd/c; the diffusion coefficients in 
the upstream and downstream region are denoted by D u and 
Dd, respectively. Fig. |21 shows the cycle-time distribution 
calculated for V u = 0.01c and r = 4, where we take D u — 
Dd = D. The solid curve represents the result obtained by 
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Figure 3. Time distribution for one-cycle of particles in a non- 
relativistic shock with V s h = 0.01c, r = 4. The solid curve is the 
approximate solution, G(i), obtained by inverting equation 1791 
numerically. The dots are results from a Monte Carlo simulation. 



inverting equation 179H numerically and the dots are the 
results from the Monte Carlo simulation. We see that the 
approximation fits the simulation results quite well. Note 
that, this distribution is very skew and peaks at t <C (At), 
which are prominent features of the diffusion process owing 
to the random walk motion of the particles. 

Thus, we can write an approximated form of (p: 



(p{Aq,s) = ifi q (Aq)G(s) 



(80) 



with the quantities in (|73^ . Obviously, the return probability 
at one-cycle is given by 



P rct = G(0) = exp(-4i/ d ) ~ 1 - 4z/ d , 



(81) 



whic h is essent ially in agreement with the previous results 
fe.g. lBelllll978ri . The first two moments of the cycle-time are 
given by 



(At) = -( — + — 



2 



< v vs v? 



(82) 



The former result coincides with the pre vious result d erived 
from the diffusion-convection equation jl)nirvll'il)S.T;. The 
state points shown in Fig.|21were also obtained by the Monte 
Carlo simulations performed above with the same parameter 
as in Fig. El 



4.2 Steady state and time-dependent solution 

Because of the isotropy of the particle distribution at the 
shock front, we have 



F 3f+ (q,t) = ±F ai (q,t), <M>+ = i- 



(83) 



Because (Aq) <C 1 and P Ict ~ 1, using the approximations 
in (1531 . we obtain 



A 



(Ag) 



(A q) ~ 4 -y^. (84) 
3 c 



Substituting the above results into equation (1491 , the steady 
state solution for large p is given by 



F&(p,oo) = 



3Q 



K - V d po 



(85) 



with cr = (r + 2)/(r - 1), where F^ ] (p)dp = F e f(q = 
\n(p/po))dq. In terms of the isotropic part of the usual phase- 
space distribution function f(p,x,t), where x is distance 
from the shock front, we have the following relation: 



/(p,0,t) 



1 



47rp : 



■F ai (p,t) = ——F s{ (q,t) 



47rp3 



We therefore obtain 
3Qo 



/(p.O.oo) 



V u — V d po 



-( ff +2) 



(86) 



(87) 



where Q' is the mono-energetic, isotropic injection rate de- 
fined for f(p, 0, t), with which (d//dt)i n j = Q'o^(p~Po) at the 
shock front. Note that, Q' has the following relation with 
Qo: 



Q'o 



The above result (I87t is equivalent to the previous results 
[e.g. equation (3.24) of |Prurvlll983l 1 ]. 

The time-dependent solution is obtained as follows. 
Adopting the approximation 1701 1. we can obtain the solu- 
tion of ^(q, s) with 



B 1 (g, 8 )~G(e a « fl ) 1 

/(<?,*) = 7^7+ ' 



B 2 (q,s)~(A q )G(e aq s) 



-dq' . 



(Aq) {Aq) J G(e a i's) 



(89) 
(90) 



Then, Q(q,s) is given by equation l)32|l . Inverting these 
results numerically, we finally obtain ^(5,4) and Q(q,t). 
Fig. 0] represents the time development of ^(q, i) calcu- 
lated for V Bh = 0.01c, r = 4, a = 1, D u = Ai = D and 
t = 500, 5000, 50000 in the unit of {At) . The self- similarity 
given by equation 15711 is evident. Fig. ^represents the time 
development of Q(q,t) calculated for the same parameters 
as in Fig. |I] except that t = 100, 500, 2500. The solutions 
given by the present method (solid curves) fit the results 
from Monte Carlo simulations (dots) quite well. 

It is interesting to make a comparison between the cut- 
off functionby our method and that by the approximation 
of iDrurvl dl99lT) . The first two cumulants ci , C2 in his ap- 
proximation [see equations (25) and (26) in his paper] are 
determined as follows. In the present situation, the diffu- 
sion coefficients are expressed as D u = = Dop/po, where 
Do is the diffusion coefficient for particles with the injection 
momentum po- From equation I82K we can determine Do in 
terms of (Af)rj as 



Do 



-(- 

4 \Vu 



(At) . 



(91) 



Using this relation, we determine for the present situation 
3c r 



ci = 



C2 



4V sh r - 1 



(p/po - 1) <At)o, 



3c 2 



r(r A + 1) 



16K 2 h (r-l)(r + l)- 



[{p/pof - 1] (At)l 



(92) 
(93) 



The factor 1/po was however missing there. 
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Figure 4. Time development of the density of state points "3/ (g, t) 
for a non-rclativistic shock with V s h = 0.01c and r = 4. We take 
a = 1, D u = Z> d and t = 500,5000,50000, where the unit of 
time is the mean cycle-time for q = 0. The self-similarity given 
by equation 1571 is evident for large q. 




Figure 5. Time development of the cut-off function @(q,t). The 
parameters are the same as in Fig. |1] except that t=100, 500, 
2500. The solid curves are numerically inverted solutions by the 
present method, and the dots are results from Monte Carlo simu- 
lations. The dotted curves are the approximate solution of Drury 
(1991). 



Subst ituting; these results into equation (22) of iDrurvl 
(1991), we obtain the approximated cut-off function repre- 
sented in Fig.[H]by dotted curves. As mentioned in his paper, 
this approximation do not provide satisfactory fits for high 
energy tail of the distributions, while it provides fairly good 
fits for lower energies [see also fig. 4 of his paper]. 



5 RELATION BETWEEN THE ONE-CYCLE 
PROPERTIES AND PARTICLE 
DISTRIBUTION 

Recently, the acceleration of particles in ultrarelativistic 
shocks has attracted some attention in relation to the 
ultrahigh-energy cosmic rays or high-energy particles at 
the external shocks associated with gamma-ray bursts. 



Many works on this issue (theories and simulations) showed 
that in such shocks the energy gain at one-cycle is ap- 
proximately given by Aq ~ In 2, owing to the charac- 
teristics of the trajectory of particles in the upstream 
region, and the power-law index is typically given by 
a ~ 2.2 (e.g. iBednarz fc OstrowskJll99a iKirk et al]l200a 
lAchterberg et alJl200lF ~ The deviation of Aq, o& q , becomes 
the same order of (Aq). Although Peacock's formula 1501 
has been often u sed in the literature to calculate the power- 
law index (e.g. iKato fc Takaharal 1200 ll : lAchterberg et ail 
2001), it can deviate from the true value because the cri- 
terion 1511 may not hold in this situation. In addition, it 
has been shown that the cycl e-time can be dominated by 
the upstream residenc e time iGallant fc Achterberd Il999l : 
lAchterberg et ai]l200lf) . Because the distribution of the up- 
stream residence time in ultrarelativistic shocks can be much 
different from that of non-relati vistic shocks (cf. Fig. [3] and 
fig. 10 of lAchterberg et alJl200ll) . the feature of the cut-off 
function can be also different from that of non-relativistic 
shocks considerably. 

To give some insights into these issues, we briefly in- 
vestigate here how the properties of the one-step probabil- 
ity density function, that is, the stochastic properties of Aq 
and At, affect the shape of the resulting particle distribu- 
tion, for several toy models. We assume the scaling law 1361 
discussed in Section [3] and the following separable form 



<p(Aq,At) = tp q (Aq)<pt{At). 



(94) 



For the following all models, we adopt P rot = 0.5, (Aq) = 
In 2, {At)o — 1 and a = 1. First, in order to investigate the 
influence of the dispersion of Aq, we consider the following 
two simple models: 



H(Aq)-H{Aq-2{Aq)) 



(model 1), and 
<p q (Aq) = — 



2{Aq) 



{Aq)/2)-H(Aq-3{Aq)/2) 



(95) 



(96) 



(model 2), where H(x) is the unit step function: H(x) = 
for x < 0, H(x) = 1 for x > and H(x) = 1/2 for x = 0. 
The standard deviation of Aq for model 2 is a half of that 
for model 1. For both models, we adopt 



<MAt) 



H(At) - H(At - 2(At}) 
' 2(At) ' 



(97) 



These models may represent typical features of the energy- 
gain and cy cle-time in ultrare l ativis tic shocks which were 
examined bv lAchterberg et alJ <l200lT) by Monte Carlo sim- 
ulations. For model 1, solving equation 1481 numerically, we 
obtain the power-law index of A ~ 0.91, while model 2 gives 
A ~ 0.98. This result explicitly shows that the power-law 
index A is generally dependent not only on (Aq) but also on 
the dispersion (JAq, although this fact is rather evident from 
equation 148 1 . Because the criterion 15H is not satisfied for 
both models, especially for model 1, Peacock's formula 115011 . 
which gives A = 1 for both models, deviates from the true 
value, [the dispersion of A q is neglected in the derivation of 
equation 1501 and that of IPeacockl il98ll) .] This fact may 
also be responsible for the small discrep ancy in the power- 
law in dices in fig. 14 of the previous paper jKato fc Takaharal 
2001) between our results, which utilized the formula 150H . 



© 2003 RAS, MNRAS 000, 1-13 



Probabilistic description of shock acceleration 11 




10 



10" 



10" 



10" 



CT 0.5 



_ I 1 1 I I 


1 1 1 1 1 1 1 1 _ 




(a) 




\ ' s -/~"~--^ steady state " 
— 1 1 1 PH 1 — H 1 — 


— 1 — 1 — 1 — 1 — 1 — \ 


V (b) " 









10 



15 



Figure 6. The distribution function of particles at the shock 
front for steady injection (a), and the cut-off function (b). The 
time-dependent solutions are calculated at t = 10 3 . The solid 
curves denote model 1 and the dotted curves denote model 2. 



Figure 7. The same as Fig. [5] except that the dotted curves 
denote model 3. Although both models have the same mean cycle- 
time, the shape of the high energy tails are considerably different 
from each other. 



and those of lEllison. Jones fc Reynolds! l)l99(l) . which were 
determined by fitting results from Monte Carlo simulations. 
Fig. |S| (a) shows the steady state and time-dependent so- 
lutions of the distribution function of particles (in an arbi- 
trary unit) and Fig. |S| (b) shows the cut-off functions, for 
both models. These solutions are obtained by the numerical 
method given in Section [2.51 The time-dependent solutions 
are calculated at t — 10 3 . It is seen that model 2 (dotted 
curve) has a sharper cut-off than model 1 (solid curve) , ow- 
ing to the smaller dispersion of Aq than that of model 1. 

Then, we investigate the influence of the dispersion of 
At on the distribution function. We define model 3, which 
has a 'diffusive' nature on the cycle-time, by 

<p t (At) = Pretff(At)/ exp(-4V/c), (98) 

where g(t) was defined in equation 1741 . Although g(t) 
was derived from the diffusion-convection equation for the 
return-time distribution function, not for the cycle-time, it 
may partly represent essential behaviour of the cycle-time 
distribution for the present situation when the motion of 
particles is mainly diffusive. In this equation, we use the nor- 
malized distribution g(At)/ exp(— 4V/c), where V — 0.01c, 
a = V{At} and D = cV{At) /A so that (At)o = 1 [cf. equa- 
tions I76I - I78I 1. and take the return probability P rct = 0.5 
independently. ip q (Aq) is taken the same as in model 1; the 
spectral index A of model 3 therefore equals to that of model 
1. Fig. |7|(a) shows the distribution function of particles and 
(b) the cut-off function, for model 1 (solid curve) and model 
3 (dotted curve) as in Fig.|S| Although both models have the 
same mean cycle-time (At), the shapes of high-energy tails 
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Figure 8. The density of state points for model 1 (solid curve), 
model 2 (dotted curve) and model 3 (dashed curve). These are 
calculated at t = 10 3 . 

are considerably different. This is attributed to the large 
dispersion of At in model 3. 

These features may be observed from the densities of 
state points shown for the three models in Fig. |g| Note that 
all models have the same mean cycle-time (At) and these 
densities are proportional to the distribution functions of 
particles at the shock front for an impulsive injection. In 
particular for the diffusive model (model 3), on account of 
its large dispersion of At, there are particles which are ac- 
celerated effectively compared with the other models while 
inefficiently accelerated particles also exist. 
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6 CONCLUSIONS 

In this paper, we have given a new description of the first- 
order Fermi acceleration in shock waves based on a stochas- 
tic, single-particle approach. The method can treat the time 
development of the particle distribution and determine the 
power-law index exactly, even for relativistic shocks. 

In Section [5] we have formulated the mechanism from 
the point of view of stochastic process. Representing the 
properties of the one-cycle of shock crossings of particles 
by a probability density function p, we have formulated the 
acceleration process as a two-dimensional Markov process 
on the q — t plane. Introducing the density of state points 
in equation 1141 . we have derived a fundamental integral 
equation I15H for describing this stochastic process. Then we 
have related the density of state points with the distribution 
function of particles at the shock front. We have also given 
numerical methods to solve the integral equation. 

In Section |!H we have investigated the case where the 
one-step probability density function satisfies the scaling 
law given by equation 1361 . Investigating the asymptotic 
behaviour of the steady state solution, we have confirmed 
that these are characterized by the power-law distribution 
in momentum space, and derived the equation 1481 which 
determines the power-law index exactly for general cases. 
We have also shown that the time-dependent solutions sat- 
isfy the similarity law given by equation 1571 except for mo- 
mentum near the injection. Then we have derived an ap- 
proximate solution which is applicable to the acceleration in 
non-relativistic shocks. 

In Section we have applied the theory to non- 
relativistic shocks as a check of our method. We have 
confirmed that the steady state solution obtained by our 
method coincides with the well-known previous results ob- 
tained by various methods and that the time-dependent so- 
lution essentially coincides with the results of Monte Carlo 
simulations. 

Finally, we have examined the relation between the 
property of the one-step probability density function and 
the steady and time-dependent solution, in Section We 
have explicitly shown that the power-law index A generally 
depends not only on the mean energy gain (Aq) but also on 
the dispersion of Aq, a& q . This can be important for rela- 
tivistic shocks. We have shown that the distribution of At 
for the diffusion process results in a quite broad high-energy 
tail in the distribution function of particles and the cut-off 
function compared with the uniform distribution such as 
19711 even for the same mean cycle-time, because of its large 
dispersion <r At • 
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APPENDIX A: NUMERICAL INVERSION OF 
LAPLACE TRANSFORMS 

Here, we briefly explain the numerical method for inver- 
sion of Lap lace transforms used in this pape r. The metho d 
is based on lAbate fc Whitd <199Eft (see also lHosono1ll984h . 
First, the inversion formula of Laplace transform is given by 



/(*) 



1 
2^i 



7+ioo 



f{s)e s ds 



M/(7 + ii/)] 



cos(yt)dy, 



(Al) 



where 7 is a real constant taken so that the vertical line 
Re(s) = 7 lies on the right of the all poles of f(s). Applying 
the trapezoidal rule to calculate the above integral with the 
step size Ay = 7r/(2t), we obtain 



/(*) 



iaoW + ^(-l) m a m (t) 

m— 1 



(A2) 



where A relates to the numerical precision ( A = r\ In 10 for 
the precision of 1Q~ V ), and 



,(t) := Re [f(X + imH)] 



(A3) 



with X := A/(2t) and H := ir/t. The number of terms to be 
summed, n, should be chosen carefully for each problem. (In 
this paper, we choose n = 80.) In order to calculate the sum 
in equation \A2l . which is an alternating series if a m does 
not change the sign, we apply Euler's transfor mation with 
van W i.jgaarden's algorithm (see chapter 5.1 of lPress et alJ 
2002). The essence of our code implemented in C++ is as 
follows: 

#include <complex> 

using namespace std; 

typedef complex<double> cplxd; 

const double PI = 3.141592653589793238; 

double 

invLaplace (cplxd Fs(const cplxdft) , 

double t, int n, double A) { 
const int nmax = 500; 
const double X(0.5*A/t), H(PI/t) ; 
if (n > nmax) n = nmax; 

// Euler summation (van Wijngaarden' s algorithm) 
double sum, sgn, wk [nmax+1] ; 
int k = 0; 
sgn = -1.0; 

sum = 0.5*(wk[0]=0.5*real(Fs(cplxd(X,0)))) ; 
for (int m=l; m<n; m++, sgn=-sgn) { 
double nxt , cur ; 

nxt = sgn*real (Fs(cplxd(X,m*H) ) ) ; 
for (int j=0; j<=k; j++) { 

cur = wk [ j ] ; wk [ j ] = nxt ; 

nxt = . 5* (cur+nxt) ; 

} 

sum += (f abs (nxt) >f abs (cur) ) ? 

nxt : 0.5*(wk[++k]=nxt) ; 

} 

return exp(0.5*A)/t * sum; 
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Here, we use the template class complex included in the 
CH — h standard library. In the arguments of the function, 
Fs (const cplxdft s) is a user-supplied, complex- valued 
function f(s) to be inverted, t is the time t for evaluation, 
n is the number of terms to be summed in 1A2I n, and A is 
A mentioned above. 
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